{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d18bc928",
   "metadata": {},
   "source": [
    "<p align=\"center\">\n",
    "    <img src=\"https://github.com/GeostatsGuy/GeostatsPy/blob/master/TCG_color_logo.png?raw=true\" width=\"220\" height=\"240\" />\n",
    "\n",
    "</p>\n",
    "\n",
    "## Data Analytics \n",
    "\n",
    "### Measures of Central Tendency in Python \n",
    "\n",
    "#### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "\n",
    "##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4104ed24",
   "metadata": {},
   "source": [
    "### Data Analytics: Measures of Central Tendency\n",
    "\n",
    "Here's a demonstration of the calculation and visualization of measures of central tendency in Python. This demonstration is part of the resources that I include for my courses in Spatial / Subsurface Data Analytics and Geostatistics at the Cockrell School of Engineering and Jackson School of Goesciences at the University of Texas at Austin.  \n",
    "\n",
    "We will cover the following statistics:\n",
    "\n",
    "#### Measures of Centrality\n",
    "* Arithmetic Average / Mean\n",
    "* Median\n",
    "* Geometric Mean\n",
    "* Harmonic Mean\n",
    "\n",
    "#### Getting Started\n",
    "\n",
    "Here's the steps to get setup in Python with the GeostatsPy package:\n",
    "\n",
    "1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). \n",
    "2. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. \n",
    "3. In the terminal type: pip install geostatspy. \n",
    "4. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. \n",
    "\n",
    "#### Importing Packages\n",
    "\n",
    "We will need some standard packages. These should have been installed with Anaconda 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4bb0eb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "plt.rc('axes', axisbelow=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f3541c0",
   "metadata": {},
   "source": [
    "#### Making Data\n",
    "\n",
    "Let's make a synthetic dataset.\n",
    "\n",
    "* You can specify the distribution type below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "id": "866c3504",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_type = 'lognormal'\n",
    "\n",
    "if dist_type == 'lognormal':\n",
    "    x = np.random.normal(loc=2.0,scale=1.0,size=10000)\n",
    "    y = np.exp(x)\n",
    "elif dist_type == 'Gaussian':\n",
    "    y = np.random.normal(loc=5000.0,scale=1000.0,size=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cdec879",
   "metadata": {},
   "source": [
    "#### Caculating Measures of Central Tendency and Plotting\n",
    "\n",
    "The following code calculates the measures of central tendency and plots them with the histogram and CDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "id": "746eb81c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Arithmetric mean =  12.35381134140441\n",
      "Geometric mean =  7.51495597840162\n",
      "Harmonic mean =  4.542792724390956\n",
      "Median =  7.466345961119204\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtUAAAINCAYAAADx3oMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABc80lEQVR4nO3deZxe893/8dcnuySIpSQklSDaiCxECIrEVty2otStJWqpFu3d1ta76vbr3rt00WorllK1U+QurVY1KLEkGmuQhJAJYg2R1STf3x/nJCaTWa7JzDXnmpnX8/GYx3Wd7znXOe/rnFk+872+55xIKSFJkiRp7XUqOoAkSZLU1llUS5IkSc1kUS1JkiQ1k0W1JEmS1EwW1ZIkSVIzWVRLkiRJzWRRLWmViLgqIr7fjNd/PCI+iIjOLZTndxHxnfz52Iioaon15uvbPSKeb6n1lUtETIqIk4rO0VZExKYRcX9ELIiIi4rOszYiIkXE1gVte2C+/S5FbF9qyyyqpVxEzI6Ifcq8jUkRsSQiBtRo2yciZpdzuy0hIsZHxPK8aP4gIl6KiN9HxDYrl0kpvZJS6p1SWl7Cuv7V2DZTSqemlL7XQvlXK1RSSg+klD7REusuSkRckL+vr9Vq/1refkFB0Yp0CvAWsF5K6Zt1LRARO0XEXRExPyLeiYhHI+KE5m64pf/xq2P9z9T4+Vue/y5ZOf3f5dqupNJYVEutbyHwnZZYUUv1CDfB5JRSb2B9YB9gMTA1IrZr6Q0V8N7aqheA42q1HZ+3V4xW7PncAng21XNns4jYBbgXuA/YGtgI+DJwQGuEa85+SCkNzf9p7Q08AJy+cjql9MOWSylpbVhUS42IiO4R8YuIeDX/+kVEdK8x/+yIeC2fd1IJH91eDBwTEVvVs70heY/2/Lxn6pAa866KiN/mvWwLgXF5D/tZEfFkRCyMiCvyj8D/kn8Efk9EbFBjHTdHxOsR8V7+MfnQpu6TlNLylNKslNJXyIqTC/J1r/bRcd4j/WKe46WIODYihgC/A3bJe9jmN/De1hiOEhH/HRFv5e/72Brtqw2TqNkbHhH3581P5Ns8unavYgn7/ZKIuDN/L4/Ud/wa28eNrSsi9o2I5/LX/hqIRg7HY0DPldvIH3vk7TUzHRQR0/L391BEDK8x79yImJXneTYiPlNj3tYRcV+e562IuDFvX2OYQM1jkO//ByPi5xHxNnBBZD9LF0bEKxExL7LhPevky28cEX+Oj3qPH4iIOv9GRcSuEfFYnumxiNh15b4l+4fi7Pw41/XJ00+Bq1NKP0kpvZUyU1NKR5W4r2ZHxJmR/by9FxE3RkSPiOgF/AXYLD7qPd4ssk8TbomIP0bE+8D4yHrKJ+frfy0ifh0R3Ro5zg2KiC9GxPSIeDci7o6ILWrMSxFxakTMyLd5SUREPq9zfkzeiogXgf+otd4T8vUuiOxn+Us15o2NiKqI+GZEvJG/lxNqzF8nIi6KiJfzffWvvO3OiDij1naerPl9J7VFFtVS474NjAFGAiOAnYDzACJif+AbZL22WwNjS1jfXOAy4P/VnhERXYH/A/4GbAKcAVwbETWHKfwn8ANgXWDlEIojgH2BbYCDyf64/zfwMbKf86/WeP1fgMH5+h8Hri0hc0P+BOxex3vpRfYPxAEppXWBXYFpKaXpwKnkvd4ppT6NvLea+gIbA5uTFU8Tau2bOqWU9sifjsi3eWOtrKXs98+RHbMNgJl5zvo0to/rXFdEbEy2P8/L3+csYLfG3h9wDR/1Vh+fT9d8f9sDVwJfIuuZvRSYGB/9cziL7Biun+f6Y0T0y+d9j2y/bAD0B35VQp6VdgZeBDbN3+OPyb5HR5L9vGwOnJ8v+02giux7dlOy7981epsjYkPgTrLvrY2AnwF3RsRGKaXxZPv6f/PjfE+t1/YEdgFuqS9wCfsK4Chgf2AQMBwYn1JaSNbb/WqN3uNX8+UPzbfZJ8+3HPg62THeBdgb+Ep9mRoTEYeS7a/DyfbfA8D1tRY7CBid5z0K+HTefnI+b3tgR+DIWq97I5+/HnAC8POI2KHG/L5k3zebAycCl8RH/8RfCIwi+9nfEDgbWAFcDXy+Rv4R+evvbPKblyqIRbXUuGOB76aU3kgpvUlWdHwhn3cU8PuU0jMppUXkPbYl+BFwcKzZSzwG6A38OKW0LKV0L/Bn4Jgay9yRUnowpbQipbQkb/tVSmleSmku2R/UR1JK/87n30b2BxOAlNKVKaUFKaWled4REbF+ibnr8irZH8y6rAC2i4h1UkqvpZSeaWRddb232r6TUlqaUrqP7I/wUfUs1xSl7PfbUkqPppSqyQqjkfWtrIR9XN+6DgSeSSndklL6EPgF8HoJ+f9I9ulHV7KC/Y+15p8CXJpSeiT/lOFqYGn+vkkp3ZxSejXf7zcCM8j+eQT4kGxIxWYppSUppUbHwtfwakrpV/n7XJLn+HpK6Z2U0gLgh3neldvpB2yRUvowH/Ne1xCO/wBmpJSuSSlVp5SuB54j+2eyMRuQ/d17rYFlGtxXuYvz/fUO2T9jIxvZ7uSU0u35/l2c94w/nOefTVa471lC/vqcCvwopTQ939c/BEbW7K0m+96en1J6BfhnjcxHAb9IKc3J38+Paq44pXRn/qlUyn/m/sbq/0R/SPb78cOU0l3AB8An8k8Zvgh8LaU0N9+XD+U/ExOBbSJicL6OLwA3ppSWNWMfSIWzqJYatxnwco3pl/O2lfPm1JhX83m98uL818B369jWnJTSilrb27yRbcyr8XxxHdO9YdVHvT+O7KP+94HZ+TIbl5K7HpsD79RuzHvujib7g/9a/pHvJxtZV2P77918vSvVPBbNUcp+r1ncLiLfp7WVuI/rW9dq3095Udno91ReKM0kK6ZmpJRqv2YL4Jv5R//zIxtyMyDfHhFxXI3hDvOB7WrkPZtsCMqjkQ2L+WJjeWqomeNjQE+yMfgrt/PXvB2yYRkzgb/lwwzOrWedtX8eYc1jVZ93yf7R69fAMg3uq1xJ3ws1rHY8ImKbyIa6vJ5/j/yQ5v0MbgH8skbed8iOWSnfv7V/h622byPigIh4OLIhOfPJ/vGrmfXtvJCvve6NyYYhzaodNv+H+Ubg83nxfQy1Pl2R2iKLaqlxr5L90Vrp43kbZD1e/WvMG0DpfgqMI/t4tOa2BsTqY0k/TjZkZKU6T8Aq0X+SfRS9D9lHtgPz9sbG7TbkM2S942tIKd2dUtqXrIh5jmzYC9T/Hhp7bxvkw0pWqnksFpIVbSv1bWRdNZWy30vVnH38GjW+h/Jxr6V+T/2BbAjFH+qYNwf4QUqpT42vniml6/PezMuA04GN8uE4T6/Mm1J6PaV0ckppM7IhEb+J7JyBlf/cNLTPax7Pt8j+wRtaI8P6KTvpjrxn/5sppS2BQ4BvRMTedbyX2j+PUOKxyj9Nmkw2XKo+9e6rxtZP6d/XvyX7eRicUlqPbOhGc34G5wBfqpV5nZTSQyW8drXvObJ9CWTnkwC3kg3j2DT/3rirxKxvkX06Ud+5B1eTfQq4N7AopTS5hHVKFc2iWlpd18hOOlr51YVsbOJ5EfGxfMzr+Xz08fpNwAmRneTWkyZc1SOlNB+4iKwncKVHyHp6zo6IrhExluxj7Rua+b5WWpfso+y3yYqhtbpiQN4bOygifkU2jryu8eGbRsSheRG8lOxj4ZU9wfOA/rF2J2f9v4joFhG7k431vDlvnwYcHhE986LvxFqvmwdsWc86W3K/N2cf3wkMjYjD8++9r1L6Pwc3AvuRfU/WdhlwakTsHJleEfEfEbEu0Ius6HsTshPTyHqqyac/GxEr/3F8N192Rf5py1yy3sbOeQ92vSdv5p8CXEY2JneTfN2bR8Sn8+cHRXZSZADvkY07XlHHqu4iGzrwnxHRJSKOBrYlG65TirPJThY8KyI2yrc9IiJWHuuG9lVj5gEblTCcal3gfeCD/NObL5eYvT6/A74VH52sun5EfLbE194EfDUi+udjoWt+QtAN6E72vVEdEQeQfY81Kj/eVwI/i+yEzc4RsUteqJMX0SvIfgfaS612waJaWt1dZL1pK78uAL4PTAGeBJ4iO/Hs+wAppb+QnTD1T7KPrh/O17O0xO39kqx4IF/fMrJi7gCynp7fAMellJ5rxnuq6Q9kH+/OBZ6tkbdUu0TEB2QFwSSyk5dGp5SeqmPZTmQncb5K9nH0nnxUPNwLPAO8HhFvNWH7r5MVdq+SjUU+tca++TmwjKywuZo1Tw68ALg6/4h8tXHYLbzf13ofp5TeAj5LdkLf22QnOz5Y4msXp5TuSSktrmPeFLIT0n5Ntv9mAuPzec+SFTaTyfbdsFrbHA08kh/3iWRjZF/M550MnJVnHQo01jN6Tr7th/NhD/cAK08GHZxPf5Bn+U1K6Z91vJe3yf6Z+ma+3bOBg/J916i893av/OvFiHgHmED2s9/gviph3c+R/RP+Yv59Vt/QpDPJPtFYQFbE31jPciVJKd0G/AS4Id+vT1P6JQIvA+4GniD73fanGutdQPaP3U1k++I/yb4HSnUm2e/Mx8h+B/yE1euOP5B9v9U+B0Bqk6Lu80AkrY3ILhf3NNC91jhDSVINEXEccEpK6VNFZ5Fagj3VUjNFxGciu/7uBmQ9Mf9nQS1J9cuHy32F7FMCqV2wqJaa70tk13KdRTaUo7njIyWp3crH0b9JNtzouoLjqAOKiCsju2HR0/XMj4i4OCJmRnZjoh3qWm6N1zn8Q5IkSR1FROxBdv7GH1JK29Ux/0Cym4AdSHYTq1+mlHZubL32VEuSJKnDSCndTx33V6jhULKCO6WUHgb6xEd3ma1Xl5YKWIRevXqljTduzvXyJWntreiUXfGt04qW7Z/oll9IblkLrnbFijxrJ/tSWsvKT4IjguXLl/Phhx+yYsWK1b569+5Nly5dWLx4Me+99x4pJVasWMHy5dlFgfr27Uu3bt1YsGAB77777hrb2GyzzejSpQvvvfce77333hrzN998czp37sz8+fN5//336dSpE/37919jOak9eeWVVxLZ1WxWmpBSasr4/c1Z/aZIVXlbQ3djbdtFNcDLL9e+sZYktY6xV40FYNL4SS263ml9svWOnN9y6x07NlvnpEktt059ZNasWVx66aVUVVVRVVXF3LlzmTt3Lrfeeiv/8R//wZ///GcOPnj1O6lHBLfddht77bUX9957Lz/60Y9Yb731WG+99dhkk01Yb731OPHEE+nbty8vvfQSzz33HOuuuy69e/ema9eudOrUia233pquXbuuVjR37tyZTp060a1bN9Zff306depESonsEuBS+xcRi1NKO7b2dtt8US1JUrktWbKEe+65hyeeeILXXnuN119/nTlz5vCVr3yF448/niVLlnDxxRez+eab079/f8aMGcPmm2/Oxz+e3aBwl112YdKkSay//vqsv/76bLDBBqy33nqrPjnYa6+92Guvverd/qBBgxg0aFC98/v06UOfPn3qnW9BLTXJXFa/02h/Srhrq0W1JEk1pJR45plnuPPOOxk2bBgHHnggS5Ys4ZBDDiGlxAYbbEDfvn3p168f666b3Whx2223ZfHixfUWrxtttBF77rlna74NSWtvInB6fqfVnYH3UkoNDv0Ai2pJkoCsmL7qqqv40Y9+xIwZM4Csh/nAAw+kT58+TJ48me22245evXqt8Vp7gqW2IyKuB8YCG0dEFfA/QFeAlNLvyO6weiDZHVUXASeUsl6LakmSgHPOOYef/vSnbLXVVlxyySXsueeefOITn1g1f+edG72illSnDz/8kKqqKpYsWVJ0lHalR48e9O/fn65duzbpdSmlYxqZn4DTmprHolqS1CEtX76cSy+9lFGjRrHzzjtz0kknsc022/DFL37Rq6SoRVVVVbHuuusycOBAP9VoISkl3n77baqqqho836A1+VtDktThTJs2jV133ZXTTjuNP/7xj6SU2GabbTjppJMsqNXilixZwkYbbWRB3YIigo022qiiev/9zSFJ6jAeeOAB9thjD7bffnumT5/Otddey8UXX2yxo7Lze6zlVdo+dfiHJKnDmDVrFrNmzeKiiy7ihBNOYIMNNig6kqR2wp5qSVK7dtNNN/HrX/8agKOOOorp06fzjW98w4JaquWkk07i2WefBeCHP/zhqvbZs2ez3XbblWWbNbcDsOuuu5ZlO63BolqS1C4tX76c//3f/+Vzn/scf/rTnwDo2bMn6623XsHJpMqzfPlyLr/8crbddltgzWK3XGpv56GHHmqV7ZaDRbUkqd15//33OeywwzjnnHM4+OCDuf3224uOJBXqsMMOY9SoUQwdOpQJEyYA0Lt3b775zW8yYsQIJk+ezNixY5kyZQrnnnsuixcvZuTIkRx77LFAVnSffPLJDB06lP3224/FixcDMHbsWL7+9a+z4447MmTIEB577DEOP/xwBg8ezHnnnbdq+3/84x/ZaaedGDlyJF/60pdYvnx5ndvp3bv3qtf85Cc/YdiwYYwYMYJzzz23tXbVWnNMtSSpXamurmavvfbiiSee4Fe/+hWnnXZaxZ3QpI5t7NixLbq+SZMmNbrMlVdeyYYbbsjixYsZPXo0RxxxBAsXLmTnnXfmoosuWm3ZH//4x/z6179m2rRpQDb8Y8aMGVx//fVcdtllHHXUUdx66618/vOfB6Bbt25MmTKFX/7ylxx66KFMnTqVDTfckK222oqvf/3rvPHGG9x44408+OCDdO3ala985Stce+21a2ynpr/85S/ccccdPPLII/Ts2ZN33nmnubup7CyqJUntSpcuXTj11FMZMGAAn/70p4uOI1WEiy++mNtuuw2AOXPmMGPGDDp37swRRxxR0usHDRrEyJEjARg1ahSzZ89eNe+QQw4BYNiwYQwdOpR+/foBsOWWWzJnzhz+9a9/MXXqVEaPHg3A4sWL2WSTTRrc3j333MMJJ5xAz549Adhwww1Lfq9FsaiWJLUL119/PZ07d+aoo47ipJNOKjqOVK9SepZbenv33HMPkydPpmfPnowdO5YlS5bQo0cPOnfuXNI6unfvvup5586dVw3/qDmvU6dOqy3XqVMnqqurSSlx/PHH86Mf/aiF3lFlcky1JKnN+8EPfsB//ud/ctVVVxUdRao47733HhtssAE9e/bkueee4+GHH270NV27duXDDz9ske3vvffe3HLLLbzxxhsAvPPOO7z88ssNbmfffffl97//PYsWLVr1mkpnUS1JarNSSlx00UWcd955HHnkkauu8iHpI/vvvz/V1dUMGTKEc889lzFjxjT6mlNOOYXhw4evOoGwObbddlu+//3vs99++zF8+HD23XdfXnvttQa3s//++3PIIYew4447MnLkSC688MJm5yi3SCkVnWGt9erVKy1cuLDoGJI6qLFXjQVg0vhJLbreaX2y9Y6c33LrXXliVGt/7Fxup59+Opdccgl77LEHf/vb31b76FmqFNOnT2fIkCFFx2iX6tq3EbEopdSrtbPYUy1JarOGDBnCWWedxT/+8Q8LakmF8kRFSVKbMm/ePB599FEOPvhgTjvttKLjSBJgUS1JakOWLFnCvvvuy5w5c5g9ezbrr79+0ZEkCbColiS1IaeffjpPPfUUd955pwW1pIrimGpJUptw7733csUVV/DVr36VAw88sOg4krQai2pJUsV74YUX2G+//RgwYADf+973io4jSWuwqJYkVbxtttmGv/71r0ybNo311luv6DiS1sIvfvGLVTdzqctJJ53Es88+24qJWlZZi+qI6BMRt0TEcxExPSJ2iYgNI+LvETEjf9wgXzYi4uKImBkRT0bEDuXMJkmqfN/97nf55S9/CcA+++zDhhtuWHAiSWuroaJ6+fLlXH755Wy77batnKrllLun+pfAX1NKnwRGANOBc4F/pJQGA//IpwEOAAbnX6cAvy1zNklSBbv33nv57ne/y4MPPlh0FKnN+973vscnPvEJPvWpT3HMMcdw4YUXMmvWLPbff39GjRrF7rvvznPPPQfA7Nmz2WuvvRg+fDh77703r7zyCgDjx4/ny1/+MmPGjGHLLbdk0qRJfPGLX2TIkCGMHz9+1bb+9re/scsuu7DDDjvw2c9+lg8++ICLL76YV199lXHjxjFu3DgAevfuzTe/+U1GjBjB5MmTGTt2LFOmTAHgr3/9KzvssAMjRoxg7733bt2dtZbKdvWPiFgf2AMYD5BSWgYsi4hDgbH5YlcDk4BzgEOBP6TsFo8P573c/VJKr5UroySpMr366qscdNBBbLPNNlx++eVFx5FaVH6D0xbT2I1SH3vsMW699VaeeOIJPvzwQ3bYYQdGjRrFKaecwu9+9zsGDx7MI488wle+8hXuvfdezjjjDI4//niOP/54rrzySr761a9y++23A/Duu+8yefJkJk6cyCGHHMKDDz7I5ZdfzujRo5k2bRr9+/fn+9//Pvfccw+9evXiJz/5CT/72c84//zz+dnPfsY///lPNt54YwAWLlzIzjvvzEUXXbRa3jfffJOTTz6Z+++/n0GDBvHOO++07A4rk3JeUm8Q8Cbw+4gYAUwFvgZsWqNQfh3YNH++OTCnxuur8rbViuqIOIWsJ5tu3bqVLbwkqRgrVqzga1/7GkuXLuXmm292DLXUTA8++CCHHnooPXr0oEePHhx88MEsWbKEhx56iM9+9rOrllu6dCkAkydP5k9/+hMAX/jCFzj77LNXLXPwwQcTEQwbNoxNN92UYcOGATB06FBmz55NVVUVzz77LLvtthsAy5YtY5dddqkzV+fOnTniiCPWaH/44YfZY489GDRoEECbGfZVzqK6C7ADcEZK6ZGI+CUfDfUAIKWUIiI1ZaUppQnABIBevXo16bWSpMr3wgsv8Ne//pXzzjuPoUOHFh1HanGN9Sy3hhUrVtCnTx+mTZvWpNd1794dgE6dOq16vnK6urqazp07s++++3L99dc3uq4ePXrQuXPnJm2/kpVzTHUVUJVSeiSfvoWsyJ4XEf0A8sc38vlzgQE1Xt8/b5MkdSCf/OQneeKJJ7jggguKjiK1C7vtthv/93//x5IlS/jggw/485//TM+ePRk0aBA333wzACklnnjiCQB23XVXbrjhBgCuvfZadt9995K3NWbMGB588EFmzpwJZEM8XnjhBQDWXXddFixYUNI67r//fl566SWANjP8o2xFdUrpdWBORHwib9obeBaYCByftx0P3JE/nwgcl18FZAzwnuOpJanjeP755/nlL3/JihUr2HLLLYmIoiNJ7cLo0aM55JBDGD58OAcccADDhg1j/fXX59prr+WKK65gxIgRDB06lDvuyEqyX/3qV/z+979n+PDhXHPNNauuwFOKj33sY1x11VUcc8wxDB8+nF122WXVCZCnnHIK+++//6oTFRtax4QJEzj88MMZMWIERx999Nq/+VYU2XmBZVp5xEjgcqAb8CJwAlkhfxPwceBl4KiU0juR/fb8NbA/sAg4IaU0paH19+rVKy1cuLBs+SWpIWOvGgvApPGTWnS90/pk6x05v+XWOzY/M2pSJXzuXIeUErvtthvPPfcczz33HJtssknRkaQWM336dIYMGVJohg8++IDevXuzaNEi9thjDyZMmMAOO7T9qxfXtW8jYlFKqVdrZynnmGpSStOAHeuYtca1UfKrfpxWzjySpMo0ceJEJk+ezK9//WsLaqkMTjnlFJ599lmWLFnC8ccf3y4K6kpT1qJakqTGrFixgrPOOottt92WU045peg4Urt03XXXFR2h3bOoliQV6pJLLmHGjBlcf/31dO3ateg4krRWyn1HRUmSGvSpT32Kr371q6tdL1eS2hp7qiVJhaiurmbx4sVsv/32bL/99kXHkaRmsadaklSI3/72t3zta19jxYoVRUeRpGazqJYktbolS5bw85//nCeeeIJOnfxTJJXb7Nmz2W677YqOsYbf/e53/OEPfyg6Rotw+IckqdVdeumlvPTSS1x66aVFR5HUiOrqarp0KU/JeOqpp5ZlvUWwe0CSWtnAAX2JiDW+Bg7oW3S0VvHBBx/wk5/8hHHjxrHvvvsWHUfqMJYvX87JJ5/M0KFD2W+//Vi8eDGXXXYZo0ePZsSIERxxxBEsWrQIgPHjx3Pqqaey8847c/bZZzN+/Hi+/OUvM2bMGLbccksmTZrEF7/4RYYMGcL48eNXbeP6669n2LBhbLfddpxzzjmr2nv37s23v/1tRowYwZgxY5g3bx4AF1xwARdeeCEAM2fOZJ999mHEiBHssMMOzJo1q/V2Tguwp1qSWtnLVfNIF67ZHmfOa/0wBTj77LN57bXXuP7664uOIhVi5d1YW0qpd3VdeenKyy67jKOOOopbb72Vww8/nJNPPhmA8847jyuuuIIzzjgDgKqqKh566CE6d+7M+PHjeffdd5k8eTITJ07kkEMO4cEHH+Tyyy9n9OjRTJs2jU022YRzzjmHqVOnssEGG7Dffvtx++23c9hhh7Fw4ULGjBnDD37wA84++2wuu+wyzjvvvNXyHXvssZx77rl85jOfYcmSJW3ufAt7qiVJreq73/0u9913H3vuuWfRUaQOZdCgQYwcORKAUaNGMXv2bJ5++ml23313hg0bxrXXXsszzzyzavnPfvazdO7cedX0wQcfTEQwbNgwNt10U4YNG0anTp0YOnQos2fP5rHHHmPs2LF87GMfo0uXLhx77LHcf//9AHTr1o2DDjpotW3XtGDBAubOnctnPvMZAHr06EHPnj3LuDdanj3VkqRW8+6777Lxxhuzxx57FB1FKkypPcstrXv37qued+7cmcWLFzN+/Hhuv/12RowYwVVXXcWkSR9l69WrV52v79Sp02rr6tSpE9XV1Q3evKlr165ExKptV1dXt8Rbqij2VEuSWsUjjzzC9ttvz6uvvlp0FEm5BQsW0K9fPz788EOuvfbaZq1rp5124r777uOtt95i+fLlXH/99SV/IrXuuuvSv39/br/9dgCWLl26anx3W2FRLUlqFT/72c9499136d27d9FRJOW+973vsfPOO7PbbrvxyU9+slnr6tevHz/+8Y8ZN24cI0aMYNSoURx66KElv/6aa67h4osvZvjw4ey66668/vrrzcrT2iKlVHSGtdarV6+0cOHComNI6qBWnmzU1I9yI6KeExUhpcS0Ptl6R85v2nobMnZsts6aH+22pscff5xRo0bxjW98g4suuqiQDFJRpk+fzpAhQ4qO0S7VtW8jYlFKqVc9Lykbe6olSWV31llnsdFGG/Gtb32r6CiSVBaeqChJKqtHHnmEe++9l5///OdsvPHGRceRpLKwqJYkldVOO+3EPffcw6677lp0FEkqG4tqSVJZRQR777130TEkqawcUy1JKpv/+q//4qyzzio6hiSVnUW1JKksPvjgA37/+98zb17HuP26pI7NolqSVBY33HAD77//Pl/60peKjiKpBY0dO5YpU6YAcOCBBzJ//vxiA1UIx1RLklpcSokJEyaw7bbbeoKi1I7dddddRUeoGPZUS5Ja3I033shjjz3G1772NSKi6DhShzd79mw++clPMn78eLbZZhuOPfZY7rnnHnbbbTcGDx7Mo48+ysKFC/niF7/ITjvtxPbbb88dd9wBwOLFi/nc5z7HkCFD+MxnPsPixYtXrXfgwIG89dZbABx22GGMGjWKoUOHMmHChFXL9O7dm29/+9uMGDGCMWPGtNshYfZUS5Ja3O677863vvUtTjzxxKKjSJUnv8NpiynxTqkzZ87k5ptv5sorr2T06NFcd911/Otf/2LixIn88Ic/ZNttt2WvvfbiyiuvZP78+ey0007ss88+XHrppfTs2ZPp06fz5JNPssMOO9S5/iuvvJINN9yQxYsXM3r0aI444gg22mgjFi5cyJgxY/jBD37A2WefzWWXXcZ5553XgjugMlhUS5Ja3Oabb84Pf/jDomNIqmHQoEEMGzYMgKFDh7L33nsTEQwbNozZs2dTVVXFxIkTufDCCwFYsmQJr7zyCvfffz9f/epXARg+fDjDhw+vc/0XX3wxt912GwBz5sxhxowZbLTRRnTr1o2DDjoIgFGjRvH3v/+93G+1EBbVkqQWk1LitNNO4/Of/7xjqaX6lNiz3NK6d+++6nmnTp1WTXfq1Inq6mo6d+7Mrbfeyic+8Ykmr3vSpEncc889TJ48mZ49ezJ27FiWLFkCQNeuXVcNA+vcuTPV1dUt8G4qj2OqJUkt5l//+he//e1vmTFjRtFRJDXRpz/9aX71q1+RUgLg3//+NwB77LEH1113HQBPP/00Tz755Bqvfe+999hggw3o2bMnzz33HA8//HDrBa8QFtWSpBZzxRVXsO6663LkkUcWHUVSE33nO9/hww8/ZPjw4QwdOpTvfOc7AHz5y1/mgw8+YMiQIZx//vmMGjVqjdfuv//+VFdXM2TIEM4991zGjBnT2vELFyv/G2mLevXqlRYuXFh0DEkd1NirxgIwafykJr0uIkgX1tF+ZjZ8YlqfbL0j5zdtvQ0Zm58YNamMHzu///779OvXj2OPPXa1M/+ljm769OkMGTKk6BjtUl37NiIWpZR6tXYWe6olSS3i+uuvZ9GiRV7xQ1KHZFEtSWoRy5Yt46ijjmKnnXYqOooktTqLaklSizjjjDO48cYbvdmLVIe2PNy2UlXaPrWoliQ12wMPPLDaXdYkfaRHjx68/fbbFVcEtmUpJd5++2169OhRdJRVvE61JKlZ3nrrLfbdd19OPfVUfvGLXxQdR6o4/fv3p6qqijfffLPoKO1Kjx496N+/f9ExVrGoliQ1y+WXX87SpUs5+eSTi44iVaSuXbsyaNCgomOozBz+IUlaa0uWLOGXv/wle+21F0OHDi06jiQVxp5qSdJa+9Of/sTrr7/O73//+6KjSFKh7KmWJK21e+65hwEDBrDPPvsUHUWSCmVRLUkVonuX7G6L89+bz/z35hMRRAQDB/QtOlq9rrzySmbMmEGXLn7wKaljs6iWpAqxtBrShdCnR/aVLsy+Xq6aV3S0Oi1ZsgSA7t27F5xEkopnUS1JarIlS5aw5ZZb8pvf/KboKJJUESyqJUlNdvvtt/Paa68xePDgoqNIUkWwqJYkNdkVV1zBFltswd577110FEmqCBbVkqQmmTt3Lv/4xz847rjj6NTJPyOSBBbVkqQmuuKKKwA47rjjCk4iSZXDayBJUoVbeam92np268SiZSvWaN+i/6bMnvN62fJ84QtfYJNNNmHrrbcu2zYkqa2xqJakCrfyUnu1xZkr6mkv7yX4Bg0axKmnnlrWbUhSW+PwD0kqk4ED+q66gUvNr7bszjvv5JZbbiGlVHQUSaoo9lRLUpm8XDWvnp7k1s/SUi644AKqq6s58sgji44iSRXFnmpJUkmee+45pkyZwvHHH190FEmqOBbVkqSSXH/99XTq1Imjjz666CiSVHEsqiVJjUopcd1117HXXnvRr1+/ouNIUsVxTLUklWjggL68XFXjyhrj8/bv9C3rJewqwRtvvEFEcMwxxxQdRZIqkkW1JJWo9omHYxdlj/dVlfcSdpVg00035fnnn2fFijWviy1JcviHJKkRKSWWLl1KRNC5c+ei40hSs0XE/hHxfETMjIhz65j/8Yj4Z0T8OyKejIgDG1unRbUkqUEPP/wwffv25aGHHio6iiQ1W0R0Bi4BDgC2BY6JiG1rLXYecFNKaXvgc8BvGluvRbUkqUHXXnstixcvZrvttis6iiS1hJ2AmSmlF1NKy4AbgENrLZOA9fLn6wOvNrbSaMt3xRowYEC65pprio4hqYOYOnUqo/p/NP1CPrx4waswatSoRpdf1V5Fg+2L564DwDqbLy5p+doenwt1/WqPCHbYYYc1ZzRg+fLlHHnkkYwcOZL/+Z//adJrJakI48aNWwY8VaNpQkppwsqJiDgS2D+ldFI+/QVg55TS6TWW6Qf8DdgA6AXsk1Ka2tB223RR3atXr7Rw4cKiY0jqICKi7hMVz6fO23bXXn5V+5k02D7tvBEAjPz+EyUtX0r72N/AfS/WnbMhjzzyCGPGjOGPf/wjxx57bJNeK0lFiIhFKaVeDcwvpaj+BlmdfFFE7AJcAWyXUqr3bG2Hf0iS6nXLLbfQpUsXDjjggKKjSFJLmQsMqDHdP2+r6UTgJoCU0mSgB7BxQyv1knqSpHodc8wxDB48mA033LDoKJLUUh4DBkfEILJi+nPAf9Za5hVgb+CqiBhCVlS/2dBKLaolSfXaYYcdmjwOW5IqWUqpOiJOB+4GOgNXppSeiYjvAlNSShOBbwKXRcTXyU5aHJ8aGT9nUS1JqtOkSZNYtGgRBx7Y6OVZJalNSSndBdxVq+38Gs+fBXZryjrLOqY6ImZHxFMRMS0ipuRtG0bE3yNiRv64Qd4eEXFxfhHuJyPCrhFJKtAPf/hDzjnnnKJjSFKb0BonKo5LKY1MKe2YT58L/COlNBj4Rz4N2QW4B+dfpwC/bYVskqQ6fPDBB9x///3svffeRUeRpDahiKt/HApcnT+/GjisRvsfUuZhoE9+jUBJUiu7++67Wbp0KYcddljRUSSpTSh3UZ2Av0XE1Ig4JW/bNKX0Wv78dWDT/PnmwJwar63K21YTEadExJSImFJdXV2u3JLUod1xxx1suOGGfOpTnyo6iiS1CeU+UfFTKaW5EbEJ8PeIeK7mzJRSiogm3YkgvyPOBMhu/tJyUSVJkN0gZsqUKRx00EF06eL57JJUirL+tkwpzc0f34iI28jutT4vIvqllF7Lh3e8kS9eyoW4JUllFhE89dRTLFiwoOgoktRmlG34R0T0ioh1Vz4H9gOeBiYCx+eLHQ/ckT+fCByXXwVkDPBejWEikqRW1LlzZ/r06VN0DElqM8o5pnpT4F8R8QTwKHBnSumvwI+BfSNiBrBPPg3ZtQJfBGYClwFfKWM2SVIdUkrsscceXHbZZUVHkaQ2pWzDP1JKLwIj6mh/m+y2j7XbE3BaufJIkhr35JNP8sADD3DccccVHUWS2pQiLqknSapQ1157LRHBwQcfXHQUSWpTLKolSUA29OPmm29m//33Z9NNN238BZKkVSyqJUkAPP3008yePZsjjjii6CiS1OZYVEuSAFhvvfU47rjj+PSnP110FElqc7yqvyQJgC222IKrr7666BiS1CbZUy1J4rXXXmPKlCmsWLGi6CiS1CZZVEtSM3Xvkt2FsPZXW3LNNdcwevRoXn311aKjSFKb5PAPSWqmpdWQLlyzPc5s/Sxr68Ybb2TUqFH079+/6CiS1CZZVEtSBxFQZw/6ZptuxKvz3ubCC+v4z0CSVBKLaknqIBL19ai/TURw7LHHtnomSWovHFMtSWLnnXemb9++RceQpDbLolqSxLXXXlt0BElq0yyqJUlsueWWRUeQpDbNolqSOrD/uqPoBJLUPlhUS1IH9c4imPBw0SkkqX2wqJakDurmJ2Dxh0WnkKT2waJakjqoe2fCZusVnUKS2geLaknqgFKCSbNgr62LTiJJ7YNFtSR1QAuWwh5bwkHbFp1EktoHi2pJ6oDW6wE3HwdHjyw6iSS1DxbVktQBzV9cdAJJal8sqiWploED+hIRa3y1FytWwFY/grP/XHQSSWo/uhQdQJIqzctV80gXrtkeZ7Z+lnJ48rXsGtXD+hWdRJLaD3uqJamDufv57NErf0hSy7GolqQO5tanYPQA2Hz9opNIUvthUS1JHUjVfJhSBQcNKTqJJLUvjqmWpA5kvR7wx2Ng9y2LTiJJ7YtFtSR1IOv1gP/coegUktT+OPxDkjqQX/8LXn6n6BSS1P7YUy1JHcgZt8OGPWGLDYtOIkntiz3VktTBjPNSepLU4iyqJakDGbQh9Fuv6BSS1P5YVEtSB5BS9riHV/2QpLKwqJakDmBJdfa4+6Bic0hSe2VRLUkdwDpds8ejRxYaQ5LaLYtqSepAencvOoEktU8W1ZLUAUx/o+gEktS+WVRLUjs3+x1444OiU0hS+2ZRLUnt3JQ5RSeQpPbPolqS2rlJs6BTFJ1Ckto3i2pJauf+MRPW71F0Cklq3yyqJakdW7QMPt4HNu5VdBJJat8sqiWpHevZDe4+BTbz1uSSVFYW1ZLUjn2wtOgEktQxWFRLUju1+EPo+//gwklFJ5Gk9s+iWpLaqX/OhIXLYFi/opNIUvtnUS1J7dSfn4V1u8OeWxadRJLaP4tqSWqnHp8LIzeDHl2LTiJJ7Z9FtSS1Q+8thsfmwNitik4iSR1Dl6IDSJJaXqeAX38GPjWo6CSS1DFYVEtSO7RuD/jyrkWnkKSOw+EfktQO/WU6zH6n6BSS1HFYVEtSO3TYVXDJg0WnkKSOw6JaktqhZcth322KTiFJHYdFtSS1Q107w24Di04hSR2HRbUktUM7fxx6dS86hSR1HBbVktSOvLc4e9x3cLE5JKmjsaiWpHZk/XWyx1N3KTaHJHU0FtWS1A5tsm7RCSSpY7GolqR25Kf/LDqBJHVMFtWS1E6sWAHfu6fl1jdwQF8iYo2vgQP6ttxGJKmd8DblktROzHwbFixtufW9XDWPdOGa7XHmvJbbiCS1E/ZUS1I78eBLRSeQpI6r7EV1RHSOiH9HxJ/z6UER8UhEzIyIGyOiW97ePZ+emc8fWO5sktSePPQy9Fmn6BSS1DG1Rk/114DpNaZ/Avw8pbQ18C5wYt5+IvBu3v7zfDlJUok+WAp7bll0CknqmMpaVEdEf+A/gMvz6QD2Am7JF7kaOCx/fmg+TT5/73x5SVIJrv883Da+6BSS1DGVu6f6F8DZwIp8eiNgfkqpOp+uAjbPn28OzAHI57+XL7+aiDglIqZExJTq6urasyWpQ7MrQpKKUbaiOiIOAt5IKU1tyfWmlCaklHZMKe3YpYsXL5EkgHP+DAddASkVnUSSOqZyVqW7AYdExIFAD2A94JdAn4jokvdG9wfm5svPBQYAVRHRBVgfeLuM+SSp3bjrOei3nj3VklSUsvVUp5S+lVLqn1IaCHwOuDeldCzwT+DIfLHjgTvy5xPzafL596Zkn4skNebND+Dp12HcVkUnkaSOq4jrVJ8DfCMiZpKNmb4ib78C2Chv/wZwbgHZJKnNmTQrexy39dq9vnsX6rxzoiSpdK0yKDmlNAmYlD9/EdipjmWWAJ9tjTyS1J78cyb07g6j+q/d65dWU8+dE5uXS5I6Es/0k6Q2blR/2LAndO1cdBJJ6rgsqiWpjTtx56ITSJKKGFMtSWohr7wL8xcXnUKSZFEtSW3Yt+6CoT/1+tSS1BQRsX9EPB8RMyOizotjRMRREfFsRDwTEdc1tk6Hf0hSG3bvzOyqH16sQ5JKExGdgUuAfcnu7v1YRExMKT1bY5nBwLeA3VJK70bEJo2t155qSWrDXl8AnxpYdApJalN2AmamlF5MKS0DbgAOrbXMycAlKaV3AVJKbzS20mjL91cZMGBAuuaaa4qOIamdmTp1ap2Xp5tatfpl615YkT0ueLXuy9nVXr7U9sVz1wFgnc0XN7j8b//0KDfddBNX/vRMBg3o2+D6X3gTFixtwZyjRq05Q5IqwLhx45YBT9VompBSmrByIiKOBPZPKZ2UT38B2DmldHqNZW4HXiC7Q3hn4IKU0l8b2m6bLqp79eqVFi5cWHQMSe1MRNR73eaa7WMXZY/3nV//dZ7Xpn3aeSMAGPn9JxpdfuNe8MYFqw//qGv5sb+B+15swZxt+G+HpPYtIhallHo1ML+UovrPwIfAUUB/4H5gWEppfn3rdfiHJLVhf/ic46klqYnmAgNqTPfP22qqAiamlD5MKb1E1ms9uKGVWlRLUht2wJCiE0hSm/MYMDgiBkVEN+BzwMRay9wOjAWIiI2BbYAXG1qpRbUktUEPv5w9LqsuNocktTUppWrgdOBuYDpwU0rpmYj4bkQcki92N/B2RDwL/BM4K6X0dkPr9ZJ6kjqsgQP68nLVvKJjrJWrp2SPne0akaQmSyndBdxVq+38Gs8T8I38qyQW1ZI6rJer5tV7Il6le/SV7NGiWpIqg7+OJamNWbgUnnyt6BSSpJosqiWpjfnXbKheUXQKSVJNFtWS1Mb8ey508be3JFUUfy1LUhtz7l7w6vmNLydJaj0W1ZLUBn2sd9EJJEk1WVRLUhty3yw4/Cqoml90EklSTRbVktSG3PoU/OU52LBn0UkkSTVZVEtSG3LvDNhzK+jZregkkqSaLKolqY2YvxiemQefGlR0EklSbRbVktRGPJLfRXGXLYrNIUlak0W1JLURK1ZkBfXoAUUnkSTV1qXoAJKk0hwwJPuSJFUee6olqQ1ICaqXF51CklQfi2pJHcLAAX2JiNW+2pLn34R1vw23P110EklSXRz+IalDeLlqHunC1dvizGKyrI3bn4Yl1TCqf9FJJEl1KamnOiKGlTuIJKl+f3kORm4GA/oUnUSSVJdSh3/8JiIejYivRMT6ZU0kSVrDo6/AuK2LTiFJqk9JRXVKaXfgWGAAMDUirouIfcuaTJK0ypJq2MuiWpIqVsknKqaUZgDnAecAewIXR8RzEXF4ucJJkjI/PQh2906KklSxSjpRMSKGAycA/wH8HTg4pfR4RGwGTAb+VL6IkqQzxxadQJLUkFJ7qn8FPA6MSCmdllJ6HCCl9CpZ77UkqYzeWVR0AklSQ0otqv8DuC6ltBggIjpFRE+AlNI15QonSR3dnPnZ47WPFxpDktSIUovqe4B1akz3zNskSWU0eXb2uOvAIlNIkhpTalHdI6X0wcqJ/HnP8kSSJK00+eXscXi/YnNIkhpWalG9MCJ2WDkREaOAxeWJJEla6eFXsseunYvNIUlqWKm3Kf8v4OaIeBUIoC9wdLlCSZJgaTU8XlV0CklSKUoqqlNKj0XEJ4FP5E3Pp5Q+LF8sSVLXTvDvr8PQC4tOIklqTKk91QCjgYH5a3aICFJKfyhLKkkSnTrBtn2LTiFJKkWpN3+5BtgKmAYsz5sTYFEtSWVy6WTYuFfRKSRJpSi1p3pHYNuUUipnGEnSR75/D3yqAm9N3r0LRMQa7Vv035TZc14vIJEkFa/UovppspMTXytjFklSbl5aRtV7sMsWcMO0otOsbmk1pDrGeceZ81o/jCRViFKL6o2BZyPiUWDpysaU0iFlSSVJHdyTK7L7ko/ZouAgkqSSlFpUX1DOEJKk1T25YhE9usDIzYpOIkkqRamX1LsvIrYABqeU7omInoC3IpCkMnk1LWP0x6FbU67RJEkqTEl3VIyIk4FbgEvzps2B28uUSZI6vIu6DeTvXyo6hSSpVKXepvw0YDfgfYCU0gxgk3KFkiRlV9mQJLUNpf7KXppSWrbyEkoR0YXsOtWSpBZ2yYevMTct484EdVy5TpJUgUrtqb4vIv4bWCci9gVuBv6vfLEkqWNKCf68/F3mU21BLUltSKlF9bnAm8BTwJeAu4DzyhVKkjqqWW/D63zI2E7rFx1FktQEpV79YwVwWf4lSSqTSbOyx9GdehcbRJLUJCUV1RHxEnWMoU4pbdniiSSpA5s0CzaiCwOje9FRJElNUOqJijvWeN4D+CywYcvHkaTmGTigLy9Xtd3bZY/qD+s+uRHhgGpJalNKHf7xdq2mX0TEVOD8lo8kSWvv5ap5pAvXbI8zWz/L2vj6HjDtb32LjiFJaqJSh3/sUGOyE1nPtVdQlaQWtmhZ0QkkSWuj1ML4ohrPq4HZwFEtnkaSOrgh/wu3p+TwD0lqY0od/jGu3EEkqSNL+anguw2CmG5BLUltTanDP77R0PyU0s9aJo4kdUxz5mePu24BTC8yiSRpbTTl6h+jgYn59MHAo8CMcoSSpI7mwdnZ464Di0whSVpbpRbV/YEdUkoLACLiAuDOlNLnyxVMkjqSu/Le6eH94Olio0iS1kKptynfFKh5TvqyvE2S1ALOHJs9dulcaAxJ0loqtaf6D8CjEXFbPn0YcHVDL4iIHsD9QPd8O7eklP4nIgYBNwAbAVOBL6SUlkVE93w7o4C3gaNTSrOb9nYkqW0asVnRCSRJzVFST3VK6QfACcC7+dcJKaUfNvKypcBeKaURwEhg/4gYA/wE+HlKaet8XSfmy58IvJu3/zxfTpLavWlz4aZpRaeQJDVHqcM/AHoC76eUfglU5T3O9UqZD/LJrvlXAvYCbsnbrybr9QY4lI96v28B9g4v1CqpA/jj43DcDUWnkCQ1R0lFdUT8D3AO8K28qSvwxxJe1zkipgFvAH8HZgHzU0rV+SJVwOb5882BOQD5/PfIhojUXucpETElIqZUV1fXni1Jbc5jc2Ckwz8kqU0rtaf6M8AhwEKAlNKrwLqNvSiltDylNJLs6iE7AZ9cu5irrXNCSmnHlNKOXbp4p3RJbdvyFTC1CkYPKDqJJKk5Si2ql6WUEtnwDSKiV1M2klKaD/wT2AXoExErq+H+wNz8+VxgQL7+LsD6ZCcsSlK79fwbsHAZ7Ni/6CTN170LRMQaXwMH9C06miSVXaldvTdFxKVkBfHJwBeByxp6QUR8DPgwpTQ/ItYB9iU7+fCfwJFkVwA5Hrgjf8nEfHpyPv/evJCXpHbrideyx1HtoKheWg3pwjXb48x5rR9GklpZo0V1frLgjWRDN94HPgGcn1L6eyMv7QdcHRGdyXrEb0op/TkingVuiIjvA/8GrsiXvwK4JiJmAu8An1ubNyRJbcnnRsKYj8PHNyg6iSSpORotqlNKKSLuSikNIzvZsCQppSeB7etof5FsfHXt9iXAZ0tdvyS1BxEwaI1TsiVJbU2pY6ofj4jRZU0iSR3M0mo44QZ48KWik0iSmqvUonpn4OGImBURT0bEUxHxZDmDSVJ79+w8uGoKzH2v6CSSpOZqcPhHRHw8pfQK8OlWyiNJHcajr2SP7eEkRUnq6BobU307sENK6eWIuDWldEQrZJKkDuGRV2CjnrClY6olqc1rbPhHzduEb1nOIJLU0fztBRi3dXayoiSpbWusqE71PJckNcOCJdCnB4zdqugkkqSW0NjwjxER8T5Zj/U6+XPy6ZRSWq+s6SSpnVq3Bzx9VtEpJEktpcGiOqXUubWCSFJH8uFy6OpvWElqN0q9pJ4kqQXtfDH85N6iU0iSWopFtSQV4N9z7amWpPbEolqSCvKpQUUnkCS1FItqSSpAl06wXd+iU0iSWopFtSQVYHg/6Nmt6BSSpJZiUS1JBRg/uugEkqSWZFEtSQU441NFJ5AktSSLaklqRa+/3/gykqS2x6JaklrR//wte0yp2BySpJZlUS1JrejemdljRLE5JEkty6JaklrJjDdh5ltFp5AklYNFtaQ2aeCAvkTEGl+V7K/PF51AklQuXYoOIElr4+WqeaQL12yPM1s/S6kemwP91oPXPFlRktode6olqZV8Zx/4zeFFp5AklYM91ZLUSgZ/LPuSJLU/9lRLUit4aDZc+ShULy86iSQpIvaPiOcjYmZEnNvAckdERIqIHRtbp0W1JLWC6x6Hr0+Ezv7WlaRCRURn4BLgAGBb4JiI2LaO5dYFvgY8Usp6/fUuSa3g/pdg+828PrUkVYCdgJkppRdTSsuAG4BD61jue8BPgCWlrDRSG76t14ABA9I111xTdAxJBZg6dSqj+tfRXkXJ7U1Ztq72F1ZkjwtebXj5DxYu5uATv8OJR+/P5z+zT6PrXzx3HQDW2Xxxi+QEeOFNWLC0+etZ6/ZRo9acIUllMG7cuGXAUzWaJqSUJqyciIgjgf1TSifl018Adk4pnV5jmR2Ab6eUjoiIScCZKaUpDW23TRfVvXr1SgsXLiw6hqQCRES9l9Qrtb0py9bVPnZR9njf+Q0vf88LsO8EuPtk2O8Tja9/2nkjABj5/SdaJCfA2N/AfS82fz1r097jXFhavWb7Fv03Zfac19ecIUnNEBGLUkq9GpjfYFEdEZ2Ae4HxKaXZpRbVXv1Dksps2qvQKWDMFkUnKcbS6vqK8HmtH0aSYC4woMZ0/7xtpXWB7YBJ+U3F+gITI+KQhgprx1RLUpmdORbe/i6s16PoJJIk4DFgcEQMiohuwOeAiStnppTeSyltnFIamFIaCDwMNFhQg0W1JLWKPusUnUCSBJBSqgZOB+4GpgM3pZSeiYjvRsQha7teh39IUhk98jL86F648GDYeuOi00iSAFJKdwF31Wo7v55lx5ayTnuqJamM/j4D7njGnmpJau8sqiWpjB58CbbrCxvXex66JKk9sKiWpDL612zYbWDRKSRJ5WZRLUll9MFSOGJ40SkkSeVmUS1JZbTX1rD7oKJTSJLKzat/SFIZ/ePUohNIklqDPdWSVAbL6rgttySp/bKolqQyuO/F7HHy7EJjSJJaiUW1JJXBfbOyx+36FptDktQ6LKolqQxW9lSv26PYHJKk1mFRLUktbPGH8OgrRaeQJLUmi2pJamGPvAzLlhedQpLUmiyqJamFbbUx/PSgolNIklqTRbUktbABfeDMsUWnkCS1JotqSWpBS6vhT0/Bu4uKTiJJak0W1ZLUgqbMgSOu/ujqH5KkjsGiWpJa0N9egE4Buw8qOokkqTVZVEtSC7pnBoweABv1KjpJ5eveBSJija+BA7xjjqS2p0vRASSpvXjtfXj4ZfjvvYtO0jYsrYZ04Zrtcea81g8jSc1kT7UktZBbnoQVCb4wqugkkqTWZk+1JLWQ03eDI4bBZusXnUSS1NrsqZakFhJhQS1JHZVFtSS1gNfeh+Ouh8erik4iSSqCRbUktYB/z4VrpsIHS4tOIkkqgkW1JLWAB16CLp1gh/5FJ5EkFcGiWpJawB1Pw64DoXf3opNIkopgUS1JLWD6G3D0iKJTSJKKYlEtSS1gWD84eGjRKSRJRbGollTRBg7oW+etrCvNk9+EAX2KTiFJKoo3f5FU0V6umlfPraxbP0tdVlQXnUCSVAnsqZakZnjzmezx+TeKzSFJKlbZiuqIGBAR/4yIZyPimYj4Wt6+YUT8PSJm5I8b5O0RERdHxMyIeDIidihXNklqKe/Oyh4Hb1xsDklSscrZU10NfDOltC0wBjgtIrYFzgX+kVIaDPwjnwY4ABicf50C/LaM2SSpRbz3SvbYyc/9JKlDK9ufgZTSaymlx/PnC4DpwObAocDV+WJXA4flzw8F/pAyDwN9IqJfufJJUnMteReWvFN0CklSJWiVvpWIGAhsDzwCbJpSei2f9Tqwaf58c2BOjZdV5W2113VKREyJiCnV1Z4hJKk4775YdAJJUqUoe1EdEb2BW4H/Sim9X3NeSikBqSnrSylNSCntmFLasUsXL14iqTjrDYAtP110CklSJShrUR0RXckK6mtTSn/Km+etHNaRP648Z34uMKDGy/vnbZJUkXptAgN2KzqFJKkSlPPqHwFcAUxPKf2sxqyJwPH58+OBO2q0H5dfBWQM8F6NYSKSVFFefQ/efBbS8qKTSJIqQTl7qncDvgDsFRHT8q8DgR8D+0bEDGCffBrgLuBFYCZwGfCVMmaTpGa5dyY8ewMsfLPoJJKkSlC2QckppX8B9d1LeO86lk/AaeXKI0kt6eGXoVM36PWxopNIkiqBV1aVpLVw9wvQZyBE56KTSJIqgUW1JDXRGwtg5ltZUS1JElhUS1KTPZpfUX+9AQ0vJ0nqOCyqJamJDvwkPHsWrLvG7akkSR2VRbUkNVGnTjBkU+jk/ackSTmLaklqgg+Wwsk3w6OvFJ1EklRJLKolqQnufxEufwTeW1J0EklSJbGolqQm+McM6N4Fdh9UdBJJUiWxqJakJnjgJdhpAPToWnQSSVIlsaiWpCZ4bA7s94miU0iSKo1FtSQ1wbB+cNSIolNIkiqNF4SSpCZ48ptFJ5AkVSJ7qiWpBCkVnUCSVMksqiWpBM/Oyx7veaHYHJKkymRRLUklePK17LHvusXmkCRVJotqSSrB/S9mj4M/VmwOSVJlsqiWpBKsLKq7e3q3JKkOFtWS1Ih3FsFzbxSdQpJUySyqJVWEgQP6EhFrfFWClODscUWnkCRVMj/IlFQRXq6aR7pwzfY4s/Wz1LZRL/jRgfDje4tOIkmqVPZUS1IDUoL7Z8H8xUUnkSRVMotqSWrAK+/Cnr+FP04tOokkqZJZVEtSAyY+mz3us02xOSRJlc2iWpIa8Nfn4BMfg09uUnQSSVIls6iWpHosXwEPzoZPDSo6iSSp0llUS1I9psyB95bAvg79kCQ1wqJakuoxegBM+S844JNFJ+lYunehzmuWDxzQt+hoklQvr1MtSfXo1AlG9S86RceztJp6rlk+r/XDSFKJ7KmWpDpUzYczbsuGgEiS1BiLakmqw5WPwiUPwUY9i04iSWoLLKolqQ5/eR7GfBwGbVR0EklSW2BRLUl1mDIHxm5VdApJUlthUS1JdaheYVEtSSqdRbUk1WHH/rCnRbUkqUReUk+S6vDYfxWdQJLUlthTLUk1LKsuOoEkqS2yqJakGi66L3t8b3GxOSRJbYtFtSTV8Nfns8f11yk2hySpbbGolqTc7Hfg/heLTiFJaossqiUpd92/i04gSWqrLKolCUgJrpkKuw8qOonq070LRMQaXwMH9C06miR5ST1JWumKo2D5CtjjN0UnUV2WVkO6cM32OHNe64eRpFosqiUJiIBdBxadQpLUVjn8Q1KHt+RDOOaP8OYHRSeRJLWGiNg/Ip6PiJkRcW4d878REc9GxJMR8Y+I2KKxdVpUS+rw7p0JN0yDKXOKTiJJKreI6AxcAhwAbAscExHb1lrs38COKaXhwC3A/za2XotqSR3e9f+GPuvAuK2LTiJJagU7ATNTSi+mlJYBNwCH1lwgpfTPlNKifPJhoH9jK42UUosnbS0DBgxI11xzTdExJLWAqVOnMqqOX1lTqyhb+9Qq+OSGS/jMly7ggD1H8/WTjmjSOl5YkT0ueLVlMy6em915Zp3NFzdrPTW98CYsWFq+fVl4+6hRa86Q1CGNGzduGfBUjaYJKaUJKyci4khg/5TSSfn0F4CdU0qn17W+iPg18HpK6fsNbbdNF9W9evVKCxcuLDqGpBYQEfVc2aG+Kz40vz3OhOuPhWOuhfu/Artv2bR1jM37MO47v2UzTjtvBAAjv/9Es9ZT09jfwH0vlm9fFt7ehv+WSWpZEbEopdSrgfklF9UR8XngdGDPlNLShrbr1T8kdWhvL4LRA7zyhyR1IHOBATWm++dtq4mIfYBvU0JBDY6pltTBnbYbPPo16OxvQ0nqKB4DBkfEoIjoBnwOmFhzgYjYHrgUOCSl9EYpK/XPiKQOzVEDktSxpJSqyYZ03A1MB25KKT0TEd+NiEPyxX4K9AZujohpETGxntWt4vAPSR3a9j+Hf389u/mLJKljSCndBdxVq+38Gs/3aeo67amW1CEtq84et9/MglqS1HwW1ZI6pHtnZo9HDC82hySpfbColtSqBg7oS0Ss8dXabnkye9x3m1bftCSpHXJMtaRW9XLVvHqvNdxalq+Aic9kz7v7W1CS1ALsqZbU4aQEVxxVdAq1lO5dqPPTj4ED+hYdTVIHYh+NpA6nS2c4eGjRKdRSllbXd6fFea0fRlKHZU+1pA6lejl87+8w662ik0iS2hN7qiV1KPe/COffDdtuWnQSSVJ7Yk+1pA7l2n9D7+5wwCeLTiJJak8sqiV1GIuWwc1PwGeHQ89uRaeRJLUnFtWSOoyJz8CCpfD5HYpOIklqb8pWVEfElRHxRkQ8XaNtw4j4e0TMyB83yNsjIi6OiJkR8WRE+CdPUot74U0YuinsuVXRSSRJ7U05e6qvAvav1XYu8I+U0mDgH/k0wAHA4PzrFOC3ZcwlqYM6fz946kzo7Gd0kqQWVrY/LSml+4F3ajUfClydP78aOKxG+x9S5mGgT0T0K1c2SR1XAXdElyR1AK3dX7NpSum1/PnrwMqLWm0OzKmxXFXetoaIOCUipkTElOrq6vIlldRufLg8ezzjtmJzSJLar8I+BE0pJSCtxesmpJR2TCnt2KWLl9mW1Ljbnsoe99q62BySpPartYvqeSuHdeSPb+Ttc4EBNZbrn7dJUrP976Ts0WtTS5LKpbWL6onA8fnz44E7arQfl18FZAzwXo1hIpK01ua+B1Orsuc9uhabRZLUfpXzknrXA5OBT0REVUScCPwY2DciZgD75NMAdwEvAjOBy4CvlCuXpI7ljqcbX0btU/cuEBGrfQ0c0LfoWJLaqbINSk4pHVPPrL3rWDYBp5Uri6SO66SdYcimsNfvik6i1ra0GtKFq7fFmfOKCSOp3fNqrZLatW5dYJwnKEqSysyiWlJZDBzQd42P3qOVLxL9q3/BeX+B1OTrDEmS1DRek05SWbxcNW+Nj94B4szW2X5K8PP7YeuNveGLJKn87KmW1C79/QV46R34wqiik0iSOgKLaknt0nf/DgP6wFEjik4iSeoILKoltTuPvAwPzoZv7JFdVk2SpHKzqJbU7mzzMdhnMJy4U9FJJEkdhX04ktqdDXrC307xBEVJUuuxp1pSu/PM6xbUkqTWZVEtqVkq4XrUK628HvX/3F3I5iVJHZjDPyQ1S9HXo67pzunZ48Hbtv62JUkdmz3VktqFlOBbd2XPj9m+2CySpI7HolpSu3Dlo/D069nzbn4GJ0lqZRbVktqFnt1gv22KTiFJ6qgsqiW1C8dsD3efUnQKSVJHZVEtqU1bVg3XPQ4fLC06iSSpI7OoltSm3fIkHHsdPDS76CRqC7p3oc5LQA4c0LfoaJLaOE/nkdSmXTMVttgguy251Jil1dRzCch5rR9GUrtiT7WkklTSTV5Wev4NuPsFOG4UdPK3mSSpQPZUSypJJd3kZaWf3Q/dOsPpuxWXQZIksKdaUhv27Dw4aSfYZN2ik0iSOjp7qiW1WQ+cBstXFJ1CkiR7qiW1Qe8s+uh5Z3+LSZIqgH+OJLU55/0le3xjQbE5JElayaJaUpvyzOtw6cPZc8dSS5IqhUW1pDblW3dB7+5Fp5AkaXUW1ZLajPtmwf89C2eNLTqJJEmrs6iW1GY8/DIM3hi+sUfRSSRJWp1FtaQ245y94MlvQs9uRSeRJGl1FtWSKt6SD+HJV7PnPboWm0WSpLpYVEuqeEf/EUb8DP45s+gkkiTVzaJaUsWb+Az86jAYt3XRSdRede8CEbHG18ABfYuOJqmN8DblkirWY69kj8fuAKftVmwWtW9LqyFduGZ7nDmv9cNIapPsqZZUsU6/LXv8+SEQUWwWSZIaYlEtqSKllF3tA+BjvYvNIklSYyyqJVWcdxZlPdOHDys6iSRJpbGolrSagQP61nnCVmt5dxEM/Slc8UirbVKSpGbzREVJq3m5al49J2yVf9srVsCJN8GbC2HbTcu/PUmSWopFtaSK8T9/g9uehgsPgl0GFp1GkqTSOfxD6qCKHuZRl+/fAyfuBN/Ys9AYkiQ1mT3VUgdV5DCP+hy8Lfz2CC+fJ0lqe+ypllSoD5d/9PyOE6Br5+KySJK0tiyqJRVm3gIY+TO47vFs2h5qSVJbZVEtqRDvLIL9JsDsd2CjXkWnkSSpeSyqJbW6J1+FPX8Dz72RDfn49CeKTiRJUvNYVEvtXCVe5WOni7OhH3ecAPtsU2gUSZJahFf/kNq5SrnKx4oVHz1/4CswcEP4WO/WzSBJUrlYVEsqu1fehc/+AcaPzqZHf7zYPJIktTSLakllNWkmHP1HWFoNG/UsOo0kSeXhmGpJZXP0NTDud9CzKzx4Ohw1suhEUtN070Kd5yQMHNC36GiSKow91ZLK5h8z4L/3hm/vDT27FZ1Garql1dRzTsK81g8jqaLZUy21E5VwlY+q+XDTtI/ukjjnO/CDAyyo1f7Ygy2pNnuqpXai6Kt8/OkpOO56WJE+aluna+tsW2pt9mBLqs2eaknN9vnr4Iir4WO9YPIZjp2WJHU89lRLarZbnoRv7gn/bz/o1b3oNJIktT6LaklN9tLbsPn60C3/DTL9LBi0UbGZJEkqksM/JJXsvllwwGWw5Y9gz99+dEKiBbUkqaOzp1pqYwYO6MvLVa17MtTtT2ePY38L6/eAb+0Fp4yBrp1bNYYkSRXLolpqY1rjKh+vvAsPzoaxW0G/9bJL5QF8f3/4+h5eIk+qz8pL7dXWs1snFi1bsUb7Fv03Zfac11sjmqQys6iWBMDc9+AX92fPB/0wuzTebw+HU3eFL+8KZ9wO396n0IhSxav/UnsrvASf1M5ZVEsFq284R309Wy1p8mzYZSCkBDv+Al5fkLWfMw6OGgHb5fex6OzZF5IkNciiWipY/cM56uvZWrvtvLUQ7n4eHngRHp2TtX3xJph+NkTABfvBLlvAiJ/BDw9cu21Iapr6hos4LERqeyyqpVbSmicYLquGNz6AqVUftZ33F7j0YejRBT41KGu75piP5n9pl1aJJqmG+oaL9Dh3nsW21MZUVFEdEfsDvwQ6A5enlH5ccCSpxbTECYYpvwX4rLeyoRqvL4BR/WHghln7f1wOM9+GF95c87VnjYXPj8p6ozt3yra744Amvw1JrcBiWyqvxmrOiOgO/AEYBbwNHJ1Smt3QOiumqI6IzsAlwL5AFfBYRExMKT1bbDKp5SythoXL4IOlsGBp9rjS4g/h949m7W8tzNo+cxV8YRQcPgyefi0bmgGwdY0f/SuPghN2yp6/vgCGbJItv+WG8PENYP/LsnlbbZx9SWq76j8Rsu5Pwer7hMwiXB1ZiTXnicC7KaWtI+JzwE+Aoxtab8UU1cBOwMyU0osAEXEDcChgUa1C1fdHKYDUxHX1OLf+edXL4bTbsucrTwyc9Ta8vyR73nc9+O+94fv3wNWfg77rZl+DNvxoHVO/3sRAktqF+sZmQ9OKcKmDKKXmPBS4IH9+C/DriIiUUr1/+qOBea0qIo4E9k8pnZRPfwHYOaV0eq3lTgFOySd3ABY3sNrOwPJmttfV1gWobmC75VBf5nKvo5TXNLZMQ/NL3eft6TiszXo8Do3naI31lLp8U/d1Y/NKORZFHIe6crTGOjwOa2qrx6Gh+R6Hln1NEccBijkW6wCP15iekFKasHKilJozIp7Ol6nKp2fly7xV71ZTShXxBRxJNqZl5fQXgF83c50TmtteT9uUAvZPnZnLvY5SXtPYMg3Nb8I+bzfHYW3W43FoO8dhbfZ1Y/NKORZFHIeWOhYeh457HJq6vz0Oa/+aIo5DkceikffaaM0JPA30rzE9C9i4ofVW0tVn5wI1T5vqn7c1x/+1QHt9y7a2lsixNuso5TWNLdPQ/FL3eXs6DmuzHo/D6ir5ODS23NrMa+/HwuPQfG31ODQ03+PQsq/pSMehMaXUnKuWiYguwPpkJyzWq5KGf3QBXgD2JnsjjwH/mVJ6ptBgdYiIKSmlHYvO0dF5HCqDx6EyeBwqg8ehMngcKkclHotSas6IOA0YllI6NT9R8fCU0lENrbdiTlRMKVVHxOnA3WRjcq6sxII6N6HxRdQKPA6VweNQGTwOlcHjUBk8DpWj4o5FfTVnRHyXbLjKROAK4JqImAm8A3yusfVWTE+1JEmS1FZV0phqSZIkqU2yqJYkSZKayaJakiRJaiaLakmSJKmZLKpbUEQcFhGXRcSNEbFf0Xk6qojYMiKuiIhbis7S0UREr4i4Ov85OLboPB2ZPweVwb8LlSEihkTE7yLiloj4ctF5OrL878SUiDio6CwtzaI6FxFXRsQb+W0pa7bvHxHPR8TMiDi3oXWklG5PKZ0MnAocXc687VULHYcXU0onljdpx9HEY3I4cEv+c3BIq4dt55pyLPw5KJ8mHgf/LpRJE4/D9JTSqcBRwG5F5G2v1uLv9jnATa2bsnVYVH/kKmD/mg0R0Rm4BDgA2BY4JiK2jYhhEfHnWl+b1Hjpefnr1HRX0XLHQS3jKko8JmR3pZqTL7a8FTN2FFdR+rFQ+VxF04+Dfxda3lU04ThExCHAncBdrRuz3buK0v9u7ws8C7zR2iFbQ8Xc/KVoKaX7I2JgreadgJkppRcBIuIG4NCU0o+ANT62iIgAfgz8JaX0eJkjt0stcRzUsppyTIAqssJ6Gv7T3uKaeCyebeV4HUZTjkNETMe/C2XR1J+H/IYeEyPiTuC6Vg3bjjXxOPQGepEV2osj4q6U0orWzFtO/tFr2OZ81OsGWcGweQPLnwHsAxwZEaeWM1gH06TjEBEbRcTvgO0j4lvlDtdB1XdM/gQcERG/Bf6viGAdUJ3Hwp+DVlffz4R/F1pXfT8PYyPi4oi4FHuqW0OdxyGl9O2U0n+R/VNzWXsqqMGe6haVUroYuLjoHB1dSultsvGLamUppYXACUXnkD8HlcK/C5UhpTQJmFRwDOVSSlcVnaEc7Klu2FxgQI3p/nmbWpfHofJ4TCqHx6IyeBwqg8ehMnTI42BR3bDHgMERMSgiugGfAyYWnKkj8jhUHo9J5fBYVAaPQ2XwOFSGDnkcLKpzEXE9MBn4RERURcSJKaVq4HTgbmA6cFNK6Zkic7Z3HofK4zGpHB6LyuBxqAweh8rgcfhIpJSKziBJkiS1afZUS5IkSc1kUS1JkiQ1k0W1JEmS1EwW1ZIkSVIzWVRLkiRJzWRRLUmSJDWTRbUk1SEilkfEtIh4OiJujoieFZDpgog4s472zSLilvz52Ij4c/78kIg4N39+WERsuxbb/EVE7NGE5QdGxOKI+HdETI+IRyNifI35B0XEd5uaQ5IqnUW1JNVtcUppZEppO2AZcGopL4qILuWNtaaU0qsppSPraJ+YUvpxPnkY0KSiOiI2AsaklO5vYqRZKaXtU0pDyO6k9l8RcUI+707g4Er4J0WSWpJFtSQ17gFg64joFRFX5r2v/46IQwEiYnxETIyIe4F/5NO3R8TfI2J2RJweEd/IX/NwRGyYv26riPhrREyNiAci4pN5+8ER8Ui+/D0RsWmNLCMiYnJEzIiIk/PlB0bE07VD5zl+HRG7AocAP81737eKiMdrLDe45nQNRwB/rbHc7Ij4Ub6OKRGxQ0TcHRGzIqLOfzpSSi8C3wC+mk8nYBJwUKk7X5LaAotqSWpA3vN8APAU8G3g3pTSTsA4siK1V77oDsCRKaU98+ntgMOB0cAPgEUppe3Jbud7XL7MBOCMlNIo4EzgN3n7v8h6iLcHbgDOrhFpOLAXsAtwfkRs1th7SCk9BEwEzsp732cB70XEyHyRE4Df1/HS3YCptdpeSSmNJPtH4yrgSGAM8P8aiPA48Mka01OA3RvLLUltSat/TClJbcQ6ETEtf/4AcAXwEHBIjXHNPYCP58//nlJ6p8br/5lSWgAsiIj3gP/L258ChkdEb2BX4OaIWPma7vljf+DGiOgHdANeqrHeO1JKi4HFEfFPYCdgGk13OXBCRHwDODpfT239gDdrtU2s8T5613iPSyOiTz3bilrTbwCN/jMgSW2JRbUk1W1x3iO7SmTV7xEppedrte8MLKz1+qU1nq+oMb2C7HdvJ2B+7W3kfgX8LKU0MSLGAhfUmJdqLVt7ulS3Av8D3AtMTSm9Xccyi8n+caip5vuo/R7r+5uyPTC9xnSPfN2S1G44/EOSSnc3cEZeXBMR26/tilJK7wMvRcRn83VFRIzIZ68PzM2fH1/rpYdGRI/8JMKxwGMlbnIBsG6N7S8hez+/pe6hH5AVwluXuP46RcRA4EKyfxRW2gZYYwy4JLVlFtWSVLrvAV2BJyPimXy6OY4FToyIJ4BngEPz9gvIhoVMBd6q9ZongX8CDwPfSym9WuK2bgDOyk9+3Cpvu5ash/lv9bzmTrLCvam2WnlJPeAm4OKUUs3CfVy+bklqNyI7EVuS1NHkY8PXTyl9p4Fl/gUclFKa30Lb3BS4LqW0d0usT5IqhUW1JHVAEXEbsBWwV0qpdm94zeV2Jhtf/mQLbXc08GFKaVpLrE+SKoVFtSRJktRMjqmWJEmSmsmiWpIkSWomi2pJkiSpmSyqJUmSpGayqJYkSZKa6f8D1slniXD33oIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('Arithmetric mean = ', str(np.average(y)))\n",
    "print('Geometric mean = ', str(stats.gmean(y)))\n",
    "print('Harmonic mean = ', str(stats.hmean(y)))\n",
    "print('Median = ', str(np.percentile(y,50)))\n",
    "\n",
    "if dist_type == 'lognormal':\n",
    "    bins=np.logspace(start=np.log10(0.01), stop=np.log10(10000), num=100)\n",
    "elif dist_type =='Gaussian':\n",
    "    bins=np.linspace(start=0, stop=10000, num=100)\n",
    "\n",
    "bins = plt.hist(y,bins=bins,color='darkorange',edgecolor='black',alpha=1.0,zorder=10)\n",
    "    \n",
    "plt.xlabel('Permeability (mD)'); plt.ylabel('Frequency'); plt.title('Log Normal Distribution and Measures of Central Tendancy')\n",
    "\n",
    "plt.vlines(np.average(y),0,np.max(bins[0])*1.1,color='black',label='arithmetic',zorder=10)\n",
    "plt.vlines(stats.gmean(y),0,np.max(bins[0])*1.1,color='blue',label='geometric',zorder=10)\n",
    "plt.vlines(stats.hmean(y),0,np.max(bins[0])*1.1,color='green',label='harmonic',zorder=10)\n",
    "plt.vlines(np.percentile(y,50),0,np.max(bins[0])*1.1,color='red',label='median',zorder=10)\n",
    "\n",
    "if dist_type == 'lognormal':\n",
    "    plt.xscale('log')\n",
    "plt.ylim([0,np.max(bins[0])*1.1]); plt.legend(loc='upper right');\n",
    "ax1 = plt.gca()\n",
    "ax2 = ax1.twinx()\n",
    "\n",
    "y_sort = np.sort(y); Fi = np.arange(1,len(y)+1) / float(len(y))\n",
    "\n",
    "ax2.plot(y_sort,Fi,c='black',ls='--',zorder=10); ax2.set_ylim([0,1]); ax2.grid(axis='y'); ax2.set_zorder(1)\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=1.5, top=1.6, wspace=0.2, hspace=0.3); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5419194",
   "metadata": {},
   "source": [
    "Note, for the lognormal distribution, the geometric mean may be eclipsed by the median.\n",
    "\n",
    "#### Comments\n",
    "\n",
    "This was a basic demonstration of calculating and plotting measures of central tendency in Python.\n",
    "\n",
    "I have other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling and many other workflows available at [Python Demos](https://github.com/GeostatsGuy/PythonNumericalDemos) and a Python package for data analytics and geostatistics at [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy). \n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "#### The Author:\n",
    "\n",
    "### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*\n",
    "\n",
    "With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. \n",
    "\n",
    "For more about Michael check out these links:\n",
    "\n",
    "#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n",
    "\n",
    "#### Want to Work Together?\n",
    "\n",
    "I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.\n",
    "\n",
    "* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! \n",
    "\n",
    "* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!\n",
    "\n",
    "* I can be reached at mpyrcz@austin.utexas.edu.\n",
    "\n",
    "I'm always happy to discuss,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "Michael Pyrcz, Ph.D., P.Eng. Associate Professor The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, The Jackson School of Geosciences, The University of Texas at Austin\n",
    "\n",
    "#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc2bddec",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
